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Improved lattice operators: the case of the topological charge density. 
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We analyze the properties of a class of improved lattice topological charge den- 
sity operators, constructed by a smearing-like procedure. By optimizing the choice 
■ of the parameters introduced in their definition, we find operators having (i) a better 

statistical behavior as estimators of the topological charge density on the lattice, i.e. 
less noisy; (ii) a multiplicative renormalization much closer to one; (iii) a large sup- 
pression of the perturbative tail (and other unphysical mixings) in the corresponding 
^ ' lattice topological susceptibility. 
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I. INTRODUCTION. 



In QCD an important role is played by topological properties. By the axial anomaly, 
matrix elements or correlation functions involving the topological charge density operator 
q{x) can be related to relevant quantities of hadronic phenomenology. We mention the 
topological susceptibility which is determinant in the explanation of the Ua{1) problem [|I|, 
and the on-shell nucleon matrix element of q{x), which can be related to the so-called spin 
content of the nucleon 0. 

Lattice techniques represent our best source of non-perturbative calculations, however 
investigating the topological properties of QCD on the lattice is a non-trivial task. In a 
lattice theory the field is defined on a discretized set and therefore the associated topological 
properties are strictly trivial. One relies on the fact that the physical continuum topological 
properties should be recovered in the continuum limit. 

From a field theoretical point of view, i.e. considering the lattice as a regulator, difficulties 
may come from unphysical divergences proportional to powers of the cut-off, which must 
be removed and therefore make the extraction of the physical signal hard. In order to 
get reliable quantitative estimates of physical quantities, one should control the unphysical 
cut-off dependent corrections even when they disappear in the continuum limit, given that 
numerical simulations are performed at finite lattice spacings, i.e. at finite values of the 
cut-off. Such corrections may be relevant, in that the typical values of the bare coupling 
where simulations are usually performed are actually not small, but ~ 1, thus few terms 
in perturbation theory are not always reliable. 

Considering a lattice version of q{x), q^{x), the classical continuum limit must be in 
general corrected by including a renormalization function. In pure QCD, where q{x) is 
renormalization group invariant, [|| 

q,ix) ^a'Zig'Mx) + 0{a') , (1) 

where Z^g^) is a finite function of the bare coupling gQ going to one in the limit (^g — 0, but 
at (jfQ ^ 1 it may be very different from one. The finite renormalization of the widely used 
lattice operator |^ 

1 ±^ 

^^(^) = "^ W-Tr [n^z^Hp,] (2) 

{Ufj^iy^x) is the product of link variables Uf^{x) around a 1 x 1 plaquette) is quite nonnegligible: 
for SU{3) Z{gl = 1) ~ 0.18 g. 

The relation of the zero-momentum correlation of two q^^ix) operators 

X. = E(9.(^)?.(0)) (3) 

with the topological susceptibility x is further complicated by an unphysical background 
term, which eventually becomes dominant in the continuum limit. (We recall that the 
definition of x requires also a prescription to define the product of operators 0.) Indeed 

xAqD = a'Z{gl)\ + M{gl) . (4) 
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Neglecting terms 0{a^), the background term MIqq) can be written in terms of mixings 
with the unity operator (so-called perturbative tail scaling as ~ a°) and with the trace of 
the energy- momentum (scaling as ~ a"^). In the case of the operator (^ and for SU{3), 
M^Qq) is already dominant at ~ 1: it is about 85% of at (^q = 1 0. As a consequence 
the uncertainty on x can be hardly made smaller than ~ 10% by using the operator (|^) and 
the heating method to evaluate Z{gl) and M{gl) 

Another problem, which has come up in some studies concerning the lattice determination 
of the on-shell proton matrix element of q{x) Jll],^, is that the lattice operator (|]) is very 
noisy, requiring very accurate statistics and therefore expensive simulations in order to get 
a reasonable uncertainty on the final result. In view of a full QCD lattice calculation the 
search for a better estimator appears a necessary step. 

We study, within the field theoretical approach, the possibility of improving the lattice 
estimator of q{x) with respect to all the problems listed above, that is we look for local 
versions of q{x) which are less noisy, have a multiplicative renormalization closer to one, 
and whose corresponding Xl not dominated by the unphysical backgroung signal M{gD 
in the region g^ ~ 1. (Any Xl defined from a local q^{x) will eventually be dominated by 
its perturbative tail in the continuum limit. For the purpose of evaluating x it would suffice 
to have a small tail at (7g ~ 1, which should be already in the scaling region.) 



II. IMPROVED TOPOLOGICAL CHARGE DENSITY OPERATORS. 

Inspired by the widely used smearing techniques, we consider the following set of opera- 
tors defined in terms of smeared links V!i'^^{x): 



q(;\x) 



1 ±4 r 

V e Tr n^*) n^'^ 



2%2 



(5) 



where H^^l is the product of smeared links V^'^''{x) around a 1 x 1 plaquette. Such smeared 
links are constructed by the following procedure: 

t>«(x) = (l-c)V;(-i)(x) + ^ E Vt'\x)V;^^~'\x + u)Vt'\x + 



v;«(x) 



^TTVi'\xyvi'\x 



1/2 



(6) 



where V^^l{x) = Vj*^(x — z/)^ V^^\x) and therefore q^j^\x) depend on the parameter c, which 



can be tuned to optimize the properties of q 
classical continuum limit, i.e. for a 0, q 



L 



X] 



x). All these operators have the correct 
a^q{x). 



Notice that the size of q^^\x) increases with increasing the integer parameter i. Neverthe- 
less g^*-* (x) can be still considered as local operators when keeping i fixed while approaching 
the continuum limit. Also, as we shall see, by optimizing the choice of the parameter c, 
a good improvement with respect to q^^\x) = q^{x) is already achieved for small values 
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of i. For SU{2) the procedure (P) keeps the smeared hnks V^^\x) belonging to the SU{2) 
group, and it is equivalent to the smearing procedures proposed in Ref. For > 3 the 
smeared links no longer belong to the SU{N) group. 

The procedure (^) may be used to improve any local operator involving link variables. 
Smearing methods to improve lattice estimators have been already widely employed in the 
study of long distance correlations, such as large Wilson loops and hadron source operators. 

One often adopts an equivalent "Schrodinger picture" of smearing, whereby lattice op- 
erators retain their original definition, while all links in the configuration undergo transfor- 
mation d^). Full consistency of this picture would then require that V^^\x) be unitary. (As 
it stands, V^''\x) is only unitary in the case of SU{2).) Projecting a matrix V onto SU{N) 
amounts to finding X G SU {N) which minimizes Tr ((Xt - V^){X - V)), or equivalently 
maximizes Tr ^^X''^^ + V^X^ . The solution is given by: 

X = iaV-' + (kV - a^l) (7) 
where a is the real root of the equation: 

n ((c^- - a')'/' + ia) = detV (8) 

i 

and di > are the eigenvalues of {V'^VY^'^. It can be verified that the lower loop results for 
Z((7q) and the perturbative tail -P(5'o) (see Sec. PTT| ) are not modified by rendering V^^'^^x) 
unitary as above. It is worth mentioning at this point that abrupt cooling leads to exactly 
the same unitary links X, for c = 1. Indeed, cooling reassigns to each link a new value, 
X^{x) in a way as to minimize the action, i.e. maximize: Tr{X^(x)V^'^\xy + X^(xyV^'^\x)) 
at c = 1, which coincides with Eq.(^. 

For > 3, instead of projecting back onto the SU{N) group we propose last step of the 
procedure which is simpler and should retain most of the advantages of the standard 
smearing procedure. 

III. PERTURBATIVE ANALYSIS. 

We have calculated Z^^^g^) to one loop for the once-smeared operator q^^^^x) with the 
Wilson action. To carry out this calculation, q^i\x) is expanded in a Taylor series in the 
gauge field A^{x), where U^{x) = exp{igoA^{x)). In Fig. |I| we show the three diagrams 
contributing to Z^^\ We find 



zi=N 

where Iq = 0.15493. At c = we recover the non-smeared results 



, (9) 



Z = 1 - 0.5362^^ ^ Q(^g4^ A^ = 2 , 

Z = 1 - 0.9084c/g + 0{g^) for A^ = 3 , (10) 
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which do not lead to a rehable estimate of Z{gQ ~ 1). As c varies, the following extreme 
values of Z are obtained: 

Z = 1 - 0.1360^^ + 0(^^) (c = 0.6495) = 2 , 

Z = 1 - 0.2472^2 ^ Qf^gi-^ ^ 0.6774) = 3 , (11) 

In both cases, Z is quite close to unity for typical values of Qq, making the one loop estimate 
more reliable. It is noteworthy that the last step in the smearing procedure (P) turns out to 
be essential to make Z approach one for c > 0. 

For q^^\x), we have also calculated the lowest perturbative contribution to the mixing 
with the unity operator Pig^), which is the dominant part of the background term M((7q) in 
the continuum limit. The corresponding diagram is shown in Fig. ^and leads to the result: 



p(c) = 0.002867 - 0.017685c + 0.0486650^ - 0.0753620^ 

+ 0.068526c^ - 0.034433c^ + 0.007445c^ . (12) 



The minimum of this everywhere-concave polynomial is p{c = 0.872) = 1.4 x 10~^. Thus, for 
all A^, the leading order of Pig^) diminishes by more than two orders of magnitude compared 
to its non-smeared value (c = 0). 

In the presence of dynamical fermions one should take into account the fact that, unlike 
pure gauge theory, the topological charge density mixes under renormalization with 
where is the singlet axial current. The nonrenormalizability property of the anomaly 
in the MS scheme means that the anomaly equation should take exactly the same form in 
terms of bare or renormalized quantities. However the renormalization of d^jl^{x) and q{x) 



is nontrivial [0]. A renormalization group analysis leads to the following relation valid for 



all matrix elements of a lattice version q^{x) of q{x) in the chiral limit ||14|| : 

{i2Nfq^) = Y{gl){R), (13) 
where 5^(fi'o) is a finite function of g^, and 

W.(a,,;(x)«„)exp£^gjd9 (14) 



is a renormalization group invariant quantity; (9^j^(x)ij_ indicates the operator d^j^{x) 
renormalized in the MS scheme, and the function jlg) is related to the anomalous dimension 
of the continuum operators q{x), d^j^{x) in the MS scheme: 'j{g) = (l/167r^)(3c^/2)A'/(yf'* + 
O (g^)- Notice that (R) is what can be naturally extracted also from experimental data. 

In perturbation theory one finds Y{gQ) = l + {zi+yi)gQ+0{gQ), where Zi is the coefficient 
of the 0{g^) term of the finite renormalization of in the pure gauge theory (cfr. Eq. (^), 
and yi turns out to be a small number: yi = —0.0486 for A^ = 3 and A"/ = 4 
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IV. NON-PERTURBATIVE ANALYSIS BY THE HEATING METHOD. 



Estimates of the multiplicative renormalizations of the operators q^^\x) and of the 
background term in the corresponding Xl obtained using the numerical heating 

method without any recourse to perturbation theory. This method relies on the idea 
that the multiplicative renormalization Zlg^) and the background term MIqq) is produced 
by short ranged fluctuations at the scale of the cut-off a. Therefore, when using a standard 
local algorithm (for example Metropolis or heat bath) to reach statistical equilibrium, the 
modes contributing to Z and M should not suffer from critical slowing down, unlike global 
quantities, such as the topological charge, which should experience a severe form of critical 
slowing 0]. 

We applied the heating method to the operators q^^\x) for 2 = 1,2 and for a number of 
values of c in the region < c < 1. We restricted our analysis to the SU{2) pure gauge 
theory, expecting no substantial differences for = 3. The measurements were performed 
at /3 = 2.6 {qq = 1.5384...), which is a typical value for the SU{2) simulations with the 
Wilson action. The local updating was performed using the heat-bath algorithm. 

An estimate of Z can be obtained by heating a configuration Co which is an approx- 
imate minimum of the lattice action and carries a definite topological charge Qifi- Such 
a configuration has been constructed by discretizing an instanton solution in the singular 
gauge 



M^) = z^f-^\ (^i^- - 5 ' ^^^^ 



where S4 = 1 and = ia^^ and exponentiating it to define link variables U^Ax) = exp iA^{x). 
Then a few cooling steps (about 5) were performed to make the configuration smoother. On 
a lattice 14^ and choosing p = 6 we obtained an instanton-like configuration carrying a 
topological charge Qla — 0.96 (all improved operators we considered gave approximately 
the same value for Co). 

One then constructs ensembles Cn of many independent configurations obtained by heat- 
ing Cq for the same number n of updating steps, averaging Q^l' = J2x ^ {^) over C„ at fixed 

n. Let us define Qi^n = {Q'i!)cu-i the average on the ensemble C„. Fluctuations of length 
/ ~ a contributing to Z should rapidly thermalize, while the topological structure of the 
initial configuration is left (approximately) unchanged for a long time. After a few heating 
steps where the short ranged modes contributing to Z get thermalized, QL,n should show a 
plateau approximately at Z Qla- The estimates of Z^'^\(3 = 2.6) from the plateaux observed 
in the heating procedure are reported in Table |, and should be compared with the value 



Z{(3 = 2.6) = 0.25(2) for the standard operator (|^) |jTO[- The plateaux formed by the ratios 
Ql\/Ql\ starting from n ~ 6 are clearly observed in Fig. |], where data for i = 1,2 and 
c = 0.8 are plotted versus n. Checks of the stability of the background topological structure 
of the initial configuration were performed at n = 8, 10, by cooling back the configurations 
(locally minimizing the action) finding Ql ~ Ql,o after few cooling steps. 

This analysis confirms the one-loop perturbative calculations, that is the improved op- 
erators we considered have a multiplicative renormalization closer to one than that of the 
initial operator q^{x). From Z{I3 = 2.6) ^ 0.25 of q^{x), we pass, by roughly optimizing 
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with respect to the parameter c, to Z^^\c = 0.8, P = 2.6) ^ 0.57 by one improving step, 
and Z(^)(c = 0.8, /5 = 2.6) ~ 0.75 by two improving steps. For larger i we expect to get 
Z^^^ closer and closer to one, as also suggested from the results of the cooling method [jl5 



On the other hand, we should not forget that increasing the number of improving steps the 
size of the operator gW(x) increases. One should find a reasonable compromise taking into 
account the size of the lattice one can afford in the simulations. 

A comparison of the above results for i = 1 with the one-loop calculation (|^) shows that 
the contribution of the higher perturbative orders is still non-negligible, but not so relevant 
as in the case of the operator without improving. 

Another important property of the improved operators we can infer from the heating 
method results is that they are much less noisy than q^{x) at fixed background. In other 
words, in the Monte Carlo determinations of the matrix elements of q{x) the contribution of 
short ranged fluctuations to the error is largely suppressed. A quantitative idea of this fact 
may come from the quantity e^'^ = AZ^^^ Z^'^\ where AZ^*) is the typical error of the data 
in the plateau during the heating procedure described above. We indeed found for c ~ 1.0 
and for an equal number of measurements 



e(o) 
e(0) 



e 



(2) 
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An estimate of the background signal M^g"^) can be obtained by measuring XL,n = 
{Ql)£„/V on ensembles of configurations £„, constructed by heating the fiat configuration for 
the same number n of updating steps Measurements were performed on a 12^ lattice. 

The plateau showed after few heating steps {n ~ 14 in this case) by the data of x|*,n should be 
placed approximately at the value of M^^^g^), since no topological activity is detected there, 
i.e. the background is still fiat (this is checked by cooling back the heated configurations), 
while the other modes contributing to M{gQ) should be already approximately thermalized 
(for a discussion of this issue see |l^). The estimates of M^'^\(3 = 2.6) from the plateaux 
observed during heating are given in Table and should be compared with the value 
M{j3 = 2.6) = 2.10(5) X 10~^ relative to the standard operator (H). In Fig. §we plot x|*,n fo^^ 
i = 1,2 and c = 1.0 as a function of the heating step n, and compare with the corresponding 
data for the standard operator. The expected plateaux are observed from n ~ 14. 

Notice the strong suppression of the background term in the improved operators. For 
c ~ 1 the reduction is about a factor 8 when performing one improving step, and about a 
factor 30 by two improving steps. For a larger number of improving steps, the suppression 
is expected to be larger. 

The suppression of the background term in Eq. (^) together with the relevant increase 
of Z should drastically change the relative weights of the contributions to Xl ^^le relevant 
region for Monte Carlo simulations. By a standard Monte Carlo simulation at /? = 2.6 on a 
16"^ lattice, we measured fo^' ^ = 1,2 and c = 0.6,0.8, 1.0. We performed 15000 sweeps 
using an overrelaxed algorithm; this sample size is already sufficient to show the better 
properties of the operators q^^\x). Data for Q^j-e given in Table |ITT[ For comparison we 
also calculated a^x by cooling ||16|. 
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For the standard operator we found Xl — 2.21(11) x 10~^, which, due to the large 
corresponding background term M = 2.10(5) x 10~^, does not allow one to determine a^x 
at this value of (3. Instead the improved operators g[*^(x) provide, using Eq. (|), reliable 
estimates of a^x having about 10% of uncertainty, which are consistent with each other and 
are also consistent with the determination from cooling: Xcooi ~ 

1.3(2) X 10~^ although the 

latter seems to be systematically lower. This fact may be explained taking into account that 
Ql = J2xQL{^)y which is used to estimate the topological charge of cooled configurations, 
underestimates the topological charge content (for the lattice size we are working with), as 
we found out explicitly when we constructed an instanton configuration on the lattice. 

The determinations of Z and M should not be subject to relevant finite size scaling effects 
(as explictly checked in Ref. ^^), since they have their origin in short ranged fluctuations. 



Thus finite size corrections to our estimates of Z and M should be negligible. Larger 
finite size effects are expected on the topological modes, as can be argued from numerical 
studies available in the literature. For this reason the measurement of Xl, which receives 
contributions also from topological modes, was performed on a larger lattice. We should 
say that we did not perform a complete analysis of the finite size corrections to Xl^ since 
our purpose was just to show the better behavior of the improved operators q^^\x) and not 
the determination of x for the SU{2) gauge model. So we limited ourselves to a numerical 
study not requiring a super-computer. 

If the improvement for SU{3) is similar to that achieved for SU{2), using the optimal 
operator for i = 2 at Qq = 1 the unphysical term in Eq. (^) is expected to become a small 
part of the total signal, allowing a precise determination of x by the field theoretical method. 

V. CONCLUSIONS. 

We have analyzed the properties of a class of improved lattice topological charge density 
operators constructed by a smearing-like procedure. Such operators look promising for the 
lattice calculation of the on-shell proton matrix element of the topological charge density 
operator in full QCD, which is related to the so-called proton spin content. Indeed their 
use should overcome the difficulty due to the large noise observed in preliminary quenched 
studies |1TT|J5[|, and they have a multiplicative renormalization much closer to one. 

Improved operators are also expected to provide a much better determination of the 
topological susceptibility by the field theoretical method in the SU{3) gauge theory, by 
strongly reducing the unphysical background term while enhancing the term containing x 
with larger values of the multiplicative renormalization. This should allow a precise and 
independent check of the alternative cooling method determinations (see e.g. Refs. [^,0), 



whose systematic errors are not completely controlled. Furthermore the improved operators 
may also open the road to a more reliable lattice investigation of the behavior of the topolog- 
ical susceptibility at the deconfinement transition, where cooling does not give satisfactory 
results fl^ . 

The smearing-like procedure (^) may be used to improve any local operator involving 
link variables, and a renormalization study would again be called for in all cases. We hope 
to return to this issue in the future. 
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FIGURES 

FIG. 1. One-loop diagrams contributing to the multiplicative renormalization of q^l\x). 

FIG. 2. Diagram contributing to the lowest order term of the perturbative tail. 

FIG. 3. QL,n/QLfi versus n in the heating procedure of an instanton configuration. Data for 
i = 1,2 at c = 0.8 are shown. For comparison the full line represents the estimate of Z for the 
standard operator (^). 

FIG. 4. Xl," versus n in the heating procedure of the flat configuration. Data for the standard 
operator and improved operators for i = 1, 2 at c = 1.0 are shown. 
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TABLES 



TABLE L We present Z^'\p = 2.6) for i = 1,2 and various values of c, as obtained by the 
heating method for a number of ~ 750 trajectories. The errors diplayed include both a statistical 
error (determined by the typical errors of data in the plateau) and a systematic error related to 
the stability of the background configuration. The numbers in this Table should be compared with 



the value Z{f3 = 2.6) = 0.25(2) for the standard operator (|) 0. 



i 


c = 0.6 


c = 0.8 


c = 1.0 


1 


0.52(2) 


0.57(2) 


0.54(2) 


2 


0.68(2) 


0.75(2) 


0.68(2) 



TABLE IL We present M(*)(/3 = 2.6) for i = 1,2 and various values of c, as obtained by the 
heating method (1000 trajectories), i.e. from the the value of Xl plateau, observed around 

n ~ 16 for all operators considered. The numbers in this Table should be compared with the value 



M{(3 = 2.6) = 2.10(5) X 10"^ for the standard operator (|). 



i 


c = 0.6 


c = 0.8 


c = 1.0 


1 


0.60(2) X 10-5 


0.37(2) X 10-5 


0.27(2) X 10-5 


2 


0.23(2) X 10-5 


0.13(2)xl0-5 


0.07(2) X 10-5 



TABLE m. We present x^-^HP = 2.6) for i = 1,2 and various values of c, as obtained by a 
standard Monte Carlo simulation on a 16'' lattice (15000 sweeps using an overrelaxed algorithm). 
We also give the corresponding values of a^x as obtained from Eq. 0. Data for a^x must be 
compared with the cooling result: Xcool = 

1.3(2) X 10-5. 





i 


c = 0.6 


c = 0.8 


c = 1.0 




1 


1.02(7) X 10-5 


0.89(7) X 10-5 


0.74(7) X 10-5 




2 


0.93(7) xlO-5 


0.92(7)xl0-5 


0.71(6) X 10-5 


a\ = {xl - M)/Z'' 


1 


1.5(3)xl0-5 


1.6(3) X 10-5 


1.6(3) X 10-5 




2 


1.5(2)xl0-5 


1.4(2) X 10-5 


1.4(2) X 10-5 
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Fig. 1 
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